1_Dynamics/7_System_identification/Identifying discrete-time models.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas\n",
    "import numpy\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will fit an ARX model of this form to a step response.\n",
    "\n",
    "$$y(k) = a_1 y(k - 1) + a_2 y(k - 2) + b_1 u(k - 1) + b_2 u(k - 2)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>yk</th>\n",
       "      <th>uk</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>k</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>-1</th>\n",
       "      <td>0.000</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.000</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.058</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.217</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.360</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.488</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.600</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0.692</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>0.772</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>0.833</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>0.888</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>0.925</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        yk  uk\n",
       "k             \n",
       "-1   0.000   0\n",
       " 0   0.000   1\n",
       " 1   0.058   1\n",
       " 2   0.217   1\n",
       " 3   0.360   1\n",
       " 4   0.488   1\n",
       " 5   0.600   1\n",
       " 6   0.692   1\n",
       " 7   0.772   1\n",
       " 8   0.833   1\n",
       " 9   0.888   1\n",
       " 10  0.925   1"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = pandas.read_csv('../../assets/data.csv', index_col='k')\n",
    "data['uk'] = 1\n",
    "data.loc[0] = [0, 1]  # input changes at t=0\n",
    "data.loc[-1] = [0, 0]  # everything was steady at t < 0\n",
    "data = data.sort_index()\n",
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = data.yk\n",
    "u = data.uk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x6215bf828>"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEGCAYAAABrQF4qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXyV5Zn/8c9FCIR9SUiAJJiwE1G2AIpjRRHFuqBVXEvVYlGrxVFnOk5t3drfjK22tp2qLSpKFaGiVZDiVkdHqwIJi8giSQxb2EJOQlhDlnP//shipIEEOMlzznO+79fLFzk5T865osnXh/u5nus25xwiIuJfrbwuQEREmpeCXkTE5xT0IiI+p6AXEfE5Bb2IiM+19uqNExISXFpamldvLyISkZYvX17knOtxPF/jWdCnpaWRnZ3t1duLiEQkM9t8vF+jpRsREZ9T0IuI+JyCXkTE5xT0IiI+p6AXEfG5RoPezGaZWaGZrTnK82ZmvzezPDNbbWYjQ1+miIicqKac0b8ATDrG8xcBA2r+mQ48ffJliYhIqDTaR++c+8jM0o5xyGTgz6563vESM+tqZr2ccztCVGPoVJbDkqeg/IDXlYiItJhQ3DCVDGyt97ig5nP/FPRmNp3qs3769OkTgrc+Tvkfwt8frK2m5d9fRMQDoQj6hhKzwd1MnHMzgZkAmZmZLb/jSVFO9Z//ng8d4lv87UVETtrDx3+SGoqumwIgtd7jFGB7CF439AK50K67Ql5Eokoogn4h8L2a7pszgNKwXJ8HKMqFhAFeVyEi0qIaXboxs7nAeCDBzAqAB4FYAOfcH4HFwLeBPOAgcHNzFXvSinJhwAVeVyEi0qKa0nVzXSPPO+COkFXUXA7tgQOFOqMXkagTPXfGBvKq/1TQi0iUiZ6gL8qt/jNhoLd1iIi0sCgK+hxo1Rq6pXldiYhIi4qeoA/kQrd0iIn1uhIRkRYVPUGv1koRiVLREfTBKijOh/j+XlciItLioiPo92yGqnJdiBWRqBQdQV/XcaOlGxGJPqEYahb+1FopIhHKOcfW4kMs3Rgga1PxCb1GlAR9TvUws/bdva5EROSYgkFHbuF+lm0MsGxTCcs2Bti19zAAXdufWNdgdAR9IE/LNiISliqqgqzdvrc62DeWkL25mD0HKwDo2TmOsenxjEnvzpj07vTv0ZGYBxt5wQZER9AX5cJADTMTEe+VVVSxcsselm0sJmtTMcs3l3CoogqA9IQOXJjRsy7YU7q1w+zkN0nyf9DXDjOL1xm9iLS8vWUVLN9UwtKaYF9dsIeKKocZDO7ZmWtGpzI6rTuj07uR2CmuWWrwf9DXDTPThVgRaX679x0ma1MxyzZW/7N+516cg9gY47TkLkz7l76MSe/GqFO606Vdy9yp7/+gr90+UGv0IhJizjkKSg7VhXrWpmLyiw4A0C42hpGndOVfJwxkdHo3RqR2o12bGE/qjIKgz9UwMxEJmcJ9Zby7dlddsO8oLQOgS7tYRqd149ox1UsxQ5O7EBsTHrcqRUHQ52iYmYiclP2HK3lnzU7eWLWNT/KKCDpI7NSWMendGZvendHp3RmY2IlWrU7+wmlz8H/Qq7VSRE5ARVWQj3J288aq7by3bidlFUFSu7fjjnP7c9mw3vRP7BiSjpiW4O+gr6qsHmY28EKvKxGRCOCcY8WWPbyxcht/+2IHxQfK6dY+limjUrl8RG9G9ukWMeFen7+DvnaYmVorReQYvtq9nwUrt/HGqu1sKT5I29atmJiRxBUjkjl7QA/atA6PtfYT5e+gV2uliBxF4b4y3vx8BwtWbWN1QSmtDM7qn8BdEwZw4dCedGzrn3j0z3fSELVWikg9DV1UPS25Cz+9eAiXDetNYufmuWHJaz4P+lwNMxOJchVVQT7O3c3rK7++qJrSrfqi6uThyfRP7Oh1ic3O30EfyNOyjUgUqr2oumDVNhat/vqi6lWjUrhiRHLEXlQ9Uf4O+qIcddyIRJGjXVS9fHgy3xoY+RdVT5R/g/5QCRzYrTN6EZ9r6KLquH4JzJgwgAtPTaJTnG6W9G/QF9V03Ki1UsR3yiqqeGftTl5bsY1/5O4m6GBocmffX1Q9Uf4N+oC2DxTxm5xd+5i7bAuvr9zGnoMVpHRrxw/H9+fyEb3pn9jJ6/LCln+DviinZpjZKV5XIiIn4VB5FYtWb2de1laWby4hNsa48NSeXD+mD2f0jQ/b+TLhxMdBn6thZiIRbO32UuYt28obq7axr6ySvj06cP+3h/CdkcnEd2zrdXkRxb9Br9ZKkYiz/3Alb36+nXnLtvB5QSltWrfi4tN6cd2YPoxOi66WyFBqUtCb2STgd0AM8Kxz7tEjnu8DzAa61hxzn3NucYhrbbqqSgh8pdZKkQjgnOOLbaXMXbaFhau2c6C8ikFJnXjw0gyuGJFM1/ZtvC4x4jUa9GYWAzwJTAQKgCwzW+icW1fvsJ8CrzjnnjazDGAxkNYM9TbNns0QrNAZvUgY21tWwYKV25i7bCvrduylXWwMl5zei+vG9mFEaledvYdQU87oxwB5zrl8ADObB0wG6ge9AzrXfNwF2B7KIo9bUU3HjVorRcJK7R2rc5dtYdHq7ZRVBMno1ZmfXz6UycN701k9782iKUGfDGyt97gAGHvEMQ8B75rZj4AOwPkNvZCZTQemA/Tp0+d4a226utZKBb1IONhzsJy/rtjGvKwt5OzaT4c2MVwxIoXrxqRyWnIXnb03s6YEfUP/BdwRj68DXnDO/drMzgReNLOhzrngN77IuZnATIDMzMwjXyN0inKhfbyGmYl4yDnHso3FzF22hcVrdlJeGWRYalce/c5pXDqsNx18NAY43DXl33QBkFrvcQr/vDQzDZgE4Jz7zMzigASgMBRFHreiXC3biHgksP8wf12xjblZW8jffYBOca25dnQq147uQ0bvzo2/gIRcU4I+CxhgZunANuBa4PojjtkCTABeMLMhQBywO5SFHpdALgyc5Nnbi0SbYNDxWX6Al5dt4d21O6mocmSe0o0fTunPxaf1ol2bGK9LjGqNBr1zrtLM7gTeobp1cpZzbq2ZPQJkO+cWAvcCz5jZ3VQv69zknGu+pZljqRtmpjN6kea2t6yCOUu2MHfZFrYUH6RLu1imnpHGtWNSGZikkQThokmLZDU98YuP+NwD9T5eB5wV2tJOUJG2DxRpbofKq5j92Sb++H9fsedgBWPTu3PvBQO58NSexMXq7D3c+O9qSO32gVqjFwm5w5VVzFu2lT98kMfufYcZP6gH904cxGkpXbwuTY7Bf0EfyNUwM5EQq6wK8tqKAn7/fh7b9hxibHp3nrphJKPT1NkWCfwX9EW50L2vhpmJhEAw6Fj0xQ6eeC+HjUUHGJbShUevPI1/6Z+g3vcI4s+g17KNyElxzvH39YX8+t0NfLlzH4OSOjFz6igmZiQp4COQv4K+qhKK82HQRV5XIhKRnHN8khfgsXc38PnWPaQndOB31w7n0tN7a+57BPNX0NcNM9MZvcjxWr65mMfe2cCS/GJ6d4njl1eexpUjU2gdE50bavuJv4K+SNsHihyvNdtK+fW7G/hgw24SOrbloUszuG5sH9q2VpukX/gr6GuHmcX397YOkQiQV7iP37yXw+IvdtKlXSz/MWkwN447hfZt/BUL4regL8rRMDORRmwJHOS37+fwxspttIuNYcaEAdxydrpGBPuYz4Je2weKHM3O0jL+539z+UvWVmJaGbec3ZfbzulH9w7awcnvfBb0Oeq4ETlCYP9hnv7wK15cspmgc1w3pg93ntefpM5xXpcmLcQ/QX+wGA4W6YxepEbpoQqe/TifWf/YyKGKKr4zMoW7JgwgtXt7r0uTFuafoA/UDjNTa6VEt4PllTz/ySZmfpRP6aEKLj69F3efP5D+iR29Lk084p+g1z6xEuXKKqp4eekWnvowj6L95UwYnMg9Fwzk1N4aOBbt/BP0gVxoFathZhJ1gkHH6yu38fi7G9hRWsa4fvH8aeogRp3SzevSJEz4J+iLcqF7uoaZSVRZs62UBxeuZfnmEoalduXXU4Yxrn+C12VJmPFX0OtCrESJ0oMVPP7uBuYs3Uy39m147KrTuXJkiubRSIP8EfQaZiZRIhh0zF++lV++vYE9B8v53plp3D1xIF3a6W+ycnT+CPq6YWY6oxf/Wl2whwcWrGXV1j2MTuvGw5eNJaN3Z6/Lkgjgj6CvG2amjhvxn5ID5Tz27gbmLttCfIe2/ObqYVwxIllz4aXJfBL0tfvEapiZ+EdV0PGXrK386p0v2VdWyc3j0vnXiQM0k0aOmz+CPpAL7RM0zEx8Y9XWPTywYA2rC0oZk96dn08eyqCenbwuSyKUP4K+KFfLNuILgf2HeeydDfwleys9Orbld9cO57JhvbVMIyfFP0E/+NteVyFywqqCjpeXbubxd3M4cLiSH5zdlxkTBtCxrT9+RcVbkf9TVDvMTKMPJEIt31zCAwvWsHb7Xsb1i+fhy05lQJKWaSR0Ij/oNcxMIlTR/sM8+taXvLq8gJ6d4/jD9SO4+LReWqaRkIv8oNc+sRJhKquCvLRkM79+L4eyiipuO6cfPzqvPx20TCPNJPJ/sopyqoeZddUwMwl/WZuK+dkba/hy5z7OHpDAQ5edSr8eGh8szSvygz6QB937QkzkfyviX4X7ynh08Zf8deU2eneJ4+kbRjJpaE8t00iLiPx0LMrRso2ErcqqILM/28xv38vhcGWQO87txx3n9qd9m8j/1ZPI0aSfNjObBPwOiAGedc492sAxVwMPAQ743Dl3fQjrbFhVJRRvhMEXN/tbiRyvJfkBHlywlg279nHOwB48dNmppCd08LosiUKNBr2ZxQBPAhOBAiDLzBY659bVO2YA8J/AWc65EjNLbK6Cv6F2mJlaKyWM7Npbxn8tXs+CVdtJ7tqOmVNHMTEjScs04pmmnNGPAfKcc/kAZjYPmAysq3fMD4AnnXMlAM65wlAX2qDaGTdqrZQwUBV0PP/JRp54L4eKoGPGhAHcfk4/2rWJ8bo0iXJNCfpkYGu9xwXA2COOGQhgZp9QvbzzkHPu7SNfyMymA9MB+vTpcyL1flPdPrEaZibe2llaxt1/WcVn+QHOG5zIg5dmcEq8lmkkPDQl6Bv6+6Zr4HUGAOOBFOBjMxvqnNvzjS9ybiYwEyAzM/PI1zh+RTkaZiaee3/9Lv5t/ueUVQR57KrTuWpUipZpJKw0JegLgNR6j1OA7Q0cs8Q5VwFsNLMNVAd/VkiqPJpAnjpuxDNlFVU8+taXvPDpJjJ6deZ/rh+hnngJS62acEwWMMDM0s2sDXAtsPCIY94AzgUwswSql3LyQ1log4pyIEHLNtLy8gr3c8VTn/LCp5v4/lnpvH7HOIW8hK1Gz+idc5VmdifwDtXr77Occ2vN7BEg2zm3sOa5C8xsHVAF/LtzLtCchVcPMwuo40ZalHOOV7K38tDCdbRrE8OsmzI5b3CS12WJHFOT+uidc4uBxUd87oF6Hzvgnpp/WkbdMDMt3UjL2FtWwU/++gWLVu/grP7xPHH1cBI7x3ldlkijIvf2PLVWSgtasaWEGXNXsqO0jB9PGsRt3+pHq1a64CqRIYKDPlfDzKTZVQUdf/y/r/jNezn07hrH/NvOZGSfbl6XJXJcIjvoNcxMmtGuvdW98Z9+FeDSYb35f1cM1cbcEpEiNyUDuVqfl2aj3njxk8gM+qoKDTOTZnG4sor/XqzeePGXyAz6Eg0zk9DLK9zPjLkrWbdjLzeflcZ9Fw2mbWvNqZHIF5lBH9D2gRI6zjnmZxfw4MK1tGsTw3M3ZjJhiHrjxT8iM+jrWit1V6ycnPq98eP6xfPENcNJUm+8+EyEBn0udOgB7dTmJifuyN74W7/Vjxj1xosPRW7Qa31eTlD93vheXdQbL/4XmUEfyFXHjZyQ+r3xl5zei//6zmnqjRffi7ygrx1mpguxcpzq98b/6qrTmaLeeIkSkRf0dbtKaelGmka98RLtIi/o61orFfTSOPXGi0Ri0BflaJiZNEq98SJfi8Cgz4P4fhpmJkel3niRb4q8tNQwMzmGTUUHmDY7i02Bg/z7hYO47Rz1xotEVtBXVUBxvlorpUGf5hVx+5wVtDKYc8tYzugb73VJImEhsoK+ZDMEK3VGL//kpSWbeWjhWtITOvDcjaPpE9/e65JEwkZkBX3tjBu1VkqNyqogP1+0jtmfbebcQT34/XUj6KQboES+IbKCvq61UsPMBEoPVnDHyyv4R14RPzg7nfsuGqL1eJEGRFbQF+VomJkAkL97P7fMzmZryUF+deXpXD061euSRMJWhAV9npZthH/kFvHDOctpHdOKObecwZj07l6XJBLWWnldwHEJ5OqO2Cj34mebuPH5ZfTq0o4Fd5ylkBdpgsg5o68bZqagj0YVVUEefnMtLy3ZwvlDEvnttSPo2DZyfnxFvBQ5vylF2j4wWu05WM4dL6/gk7wAt57Tlx9fOFgXXUWOQwQFfW1rpTpuokle4X5umZ3F9j1lPD5lGFeNSvG6JJGIEzlBH8iFmDYaZhZFPsrZzR0vr6Bt61bMnT6WUadoPV7kRERO0BflQfe+GmYWBZxzvPDpJn6+aB0Dkzrx7I2ZpHTTna4iJypyUrMoB3oM8roKaWYVVUEeWLCWucu2cEFGEk9cM5wOuugqclIi4zeoqgJKNsKQS72uRJpRyYFybp+znCX5xfxwfD/+7YJBtNJFV5GT1qQ+ejObZGYbzCzPzO47xnFXmZkzs8zQlQiUbKoZZqbWSr/K3bWPyU9+woote/jtNcP58aTBCnmREGn0jN7MYoAngYlAAZBlZgudc+uOOK4TMANYGvIq1Vrpax9sKGTGyytpGxvDvOlnMLKPRlyIhFJTzujHAHnOuXznXDkwD5jcwHE/B34FlIWwvmpqrfQl5xzPfpzPtBeySO3engV3nqWQF2kGTQn6ZGBrvccFNZ+rY2YjgFTn3KJjvZCZTTezbDPL3r17d9OrDOTWDDPr2vSvkbBWXhnkvte+4Bd/W8/EjCRevf1Mkru287osEV9qysXYhhZKXd2TZq2AJ4CbGnsh59xMYCZAZmama+TwrxXladnGRwL7D3P7SytYtqmYH53Xn7vPH6j1eJFm1JSgLwDqz4BNAbbXe9wJGAp8aGYAPYGFZnaZcy47JFUW5ajjxic27NzHtNlZFO47zO+uHc7k4cmNf5GInJSmBH0WMMDM0oFtwLXA9bVPOudKgYTax2b2IfBvIQv5AwE4VKwzeh94f/0uZsxdSfu2rXnl1jMZnqqlOJGW0GjQO+cqzexO4B0gBpjlnFtrZo8A2c65hc1aYd2uUmqtjFTOOZ75OJ//futLTu3dmWe+l0mvLlqPF2kpTbphyjm3GFh8xOceOMqx40++rHqKFPSR7HBlFfe/voZXlxdw8Wm9eHzKMNq1ifG6LJGoEv53xmqYWcQq2n+Y215cTvbmEu6aMIC7JgzQRVcRD4R/0BflVg8za6WzwEiyqegA331uKbv3HeYP14/gktN7e12SSNSKjKBPHOx1FXIc1m3fy/dmLSPoHK/ceibDdNFVxFPhvWds7TAzbQgeMbI3FXPNzM+IjTGFvEiYCO8z+rphZmqtjAQfbCjk9peW07tLO168ZazudBUJE+Ed9LUzbtRxE/YWfr6de/6yikE9OzH7+2NI6NjW65JEpEaYB31Na6WGmYW1l5Zs5mcL1jA6rTvP3phJ57hYr0sSkXrCO+gDudAhUcPMwpRzjqc+/IrH3tnAhMGJPHnDSOJi1R0lEm7CO+iLcrVsE6acc/zX4vU88/FGLh/em8emDCM2Jryv7YtEq/D+zVTQh6XKqiD/8dpqnvl4IzeeeQq/uXq4Ql4kjIXvGX3tMDO1VoaVw5VV3DV3FW+v3cmMCQO4+/wB1EwtFZEwFb5BH9D2geFm/+FKbn0xm0/yAjxwSQbf/5d0r0sSkSYI36CvG2amjptwUHKgnJteyGLNtlJ+PWUYV45K8bokEWmiMA76HA0zCxM7S8uY+txSNhcf5I/fHcXEjCSvSxKR4xC+QR/Ig+79NMzMY7XDyUoOlPPCzaMZ1y+h8S8SkbASvkFflAOJQ7yuIqrVDierCgaZO/0MTk/R/QwikSg8e+KqKqrn3OhCrGfqDyebf9uZCnmRCBaeZ/S1w8zUWumJDzcUcttLy+nVpR0vThtDSrf2XpckIichPINew8w88+bn27n7L6sYmNSJP0/TcDIRPwjToNcwMy/MWbqZn76xhtGndOfZmzScTMQvwjfoNcysxdQfTnbe4ESevH6kNvAW8ZHwDPpAri7EthDnHP/91pfM/CifycN787iGk4n4Tnj+Rhfl6I7YFlA7nGzmR/l878xTeELDyUR8KfzO6A8E4FCJzuib2TeGk53Xn7snDtRwMhGfCr+grx1mptbKZnPgcCXTa4aT/eySDKZpOJmIr4Vf0Ne1VmrppjmUHCjn5hey+GJbKY9PGcZVGk4m4nthGPS5GmbWTOoPJ3v6hpFccGpPr0sSkRYQnkGvYWYhp+FkItEr/II+kAuJGV5X4Svrd+xl6nPVw8le/sEZDEvV/Qki0SS8eunqhpnpQmyorNq6h2v+9BmtW1UPJ1PIi0SfJgW9mU0ysw1mlmdm9zXw/D1mts7MVpvZ+2Z2YgvsxRs1zCyE1mwrZepzS+navg2v3n4m/RM7eV2SiHig0aA3sxjgSeAiIAO4zsyOXFtZCWQ6504HXgV+dULVaJ/YkPly516mPreUznGxvPyDsZpAKRLFmnJGPwbIc87lO+fKgXnA5PoHOOc+cM4drHm4BDixnj21VoZEXuF+vvvsUtq0bqWQF5EmBX0ysLXe44Kazx3NNOCthp4ws+lmlm1m2bt37/7nA4ryoGMSxHVpQlnSkE1FB7j+mSWA8fIPzuCU+A5elyQiHmtK0Dd0X7xr8ECz7wKZwGMNPe+cm+mcy3TOZfbo0eOfDyjK0fr8SSgoOcgNzy6loirInFvG0q9HR69LEpEw0JSgLwBS6z1OAbYfeZCZnQ/cD1zmnDt8QtUEctVxc4J2lB7i+meWsq+sghenjWVQT114FZFqTQn6LGCAmaWbWRvgWmBh/QPMbATwJ6pDvvCEKqkbZqagP16F+8q44ZmlFB8o58/TxjI0WUtfIvK1RoPeOVcJ3Am8A6wHXnHOrTWzR8zssprDHgM6AvPNbJWZLTzKyx1d7YVYLd0cl+ID5Xz32aXsKC3j+ZtHM1x98iJyhCbdGeucWwwsPuJzD9T7+PyTrqSutVJB31SlByv47rNL2Rw4yPM3jWZ0WnevSxKRMBQ+d8YW5UBMW+jax+tKIsK+sgq+N2speYX7+dPUUYzrr9k1ItKwMAr6PIjXMLOmOHC4kpufz2Lt9r08ecNIxg9K9LokEQlj4RP0gVyI141SjTlUXsW02Vms2FLC768bwcSMJK9LEpEwFx5BX1lePedG6/PHVFZRxfQXs1m6sZjfXD2cb5/Wy+uSRCQChEfQl2wCV6UZN8dQXhnkzpdX8HFuEb/8zulcPuJYNyeLiHwtPIJerZXHVFkV5K55K/n7+kJ+fvlQrh6d2vgXiYjUCI+NR+paK7VGf6SqoOPe+Z/z1pqd/PTiIUw9Q1ssijRVRUUFBQUFlJWVeV3KcYuLiyMlJYXY2NiTfq3wCPqiXA0za0Aw6LjvtdUsWLWdH08axC1n9/W6JJGIUlBQQKdOnUhLS8OsobFd4ck5RyAQoKCggPT09JN+vTBZusnV+vwRnHM8sHAN85cXMGPCAH44Xn/bETleZWVlxMfHR1TIA5gZ8fHxIfubiPdB71zN1EoFWS3nHD9ftJ6Xlmzh1nP6cvf5unYhcqIiLeRrhbJu74P+YADK9qi1soZzjl+9s4FZn2zkpnFp3DdpcMT+oIpIePA+6Iu0fWB9v38/j6c//Irrx/bhwUszFPIiPvPhhx9yySWXtOh7hkHQ17ZWaunm6Q+/4om/53DVqBR+MXmoQl5EQsL7rptAroaZAbP+sZFfvv0llw7rzS+vPJ1WrRTyIqH08JtrWbd9b0hfM6N3Zx689NSjPv+zn/2MhIQE7rrrLgDuv/9+kpK+HluSlZXF9OnTee211+jbt/m66sLgjF7DzF5asplHFq1j0qk9+c3Vw4hRyIv4wrRp05g9ezYAwWCQefPmkZxcfVf7p59+ym233caCBQuaNeQhHM7oi3Kg51Cvq/DMK9lb+ekbazhvcCK/v24EsTHe/79XxI+OdebdXNLS0oiPj2flypXs2rWLESNGEB8fz/r165k+fTrvvvsuvXv3bvY6vA36yvLqOTenXuFpGV5ZsGob//Haas4ekMBTN4ykTWuFvIjf3HLLLbzwwgvs3LmT73//+wD06tWLsrIyVq5c2SJB722ylGysGWYWfa2Vb6/ZwT2vfM6YtO7MnJpJXGz0Ll2J+NkVV1zB22+/TVZWFhdeeCEAXbt25W9/+xs/+clP+PDDD5u9Bm/P6Iuic/vA99fv4kdzVzI8tSuzbhpNuzYKeRG/atOmDeeeey5du3YlJubr3/WkpCTefPNNLrroImbNmsXYsWObrQaPgz76plZ+lLOb219awZBenXn+5tF0aOv9ZRIRaT7BYJAlS5Ywf/58AMaPH8/48eMB6NOnD2vXrm32GrxdugnkQceeENfZ0zJaymdfBZj+Yjb9Ejvy5++PoXPcyU+lE5HwtW7dOvr378+ECRMYMMC7E1rvl26iZNnm07wibvlzNqnd2vPStDF0bd/G65JEpJllZGSQn5/vdRleB32Orztuig+Us2DVNuZnF7Bux17SEzow55axxHds63VpIhJFvAv6YGXNMDN/zbipqAryfxt28+ryAt7/chcVVY7TU7rwyORTuXxEspZrRKTFeRf0lTVzln2ydLNh5z5eXb6V11dup2j/YRI6tuGmcWlcOSqFwT2j4xqEiIQnD4P+cPWfERz0ew6Ws/Dz7by6vIDVBaW0bmVMGJLIlFGpnDOoh+5yFZGw4O0ZfUxb6BJZG11XVgX5OLeIV5cX8N66XZRXBcno1ZkHLslg8vDeWn8XkSYZP348jz/+OJmZmc3+Xt6e0ccPjiZjMRoAAAdoSURBVJhhZnmF+5i/vIDXV2yjcN9hundoww1n9OGqUSmc2lt73YpI+PL2jD4hvGfQlx6qYNHq7czPLmDV1j3EtDLOHZTIVaNSOG9wombTiESSt+6DnV+E9jV7ngYXPXrMQzZt2sQll1zCmjVrAHj88cfZv39/3fPBYJCbb76Z1NRUfvGLX4S2vhoeBn15WN4RWxV0fJJXvTTzztqdHK4MMiipEz+9eAiThyfTo5OWZkQkNCorK7nhhhsYOnQo999/f7O9j4d99C6sWivzd+/ntRUF/HXFNnaUltGlXSzXjE5lyqhUhiZ31m5PIpGukTNvL9x6661cffXVzRry0MSgN7NJwO+AGOBZ59yjRzzfFvgzMAoIANc45zY1+sIeL93sK6vgb6t38OryArI3l9DK4JyBPfjZJRlMGJJI29aRcf1ARMJX69atCQaDdY/LysrqPh43bhwffPAB9957L3Fxcc1XQ2MHmFkM8CQwESgAssxsoXNuXb3DpgElzrn+ZnYt8Evgmkbf3YOlm2DQsSQ/wPzlBby1ZgdlFUH6J3bkvosG850RySR2br5/2SISfZKSkigsLCQQCNCxY0cWLVrEpEmTgOodqD766COmTJnC66+/TuvWzbPI0pRXHQPkOefyAcxsHjAZqB/0k4GHaj5+FfiDmZlzzh3tRStpzcSnVp5Q0Sdjz6EKdu87TKe41lw5MoUpmakMS+mipRkRaRaxsbE88MADjB07lvT0dAYPHvyN5++55x5KS0uZOnUqc+bMoVWr0Dd5NCXok4Gt9R4XAEcOTq47xjlXaWalQDxQVP8gM5sOTAfo26sb45I6nmDZJ65NTCvOG5LEBRlJ2uxDRFrEjBkzmDFjxlGff/jhh5v1/ZsS9A2d6h55pt6UY3DOzQRmAmRmZrqnbhjVhLcXEZGT0ZS/IxQA9W9fTQG2H+0YM2sNdAGKQ1GgiIicnKYEfRYwwMzSzawNcC2w8IhjFgI31nx8FfC/x1qfFxFpKZEaRaGsu9Ggd85VAncC7wDrgVecc2vN7BEzu6zmsOeAeDPLA+4B7gtZhSIiJyguLo5AIBBxYe+cIxAIhKzl0rz6F5CZmemys7M9eW8RiQ4VFRUUFBR8o3c9UsTFxZGSkkJs7Df3sDCz5c6545qEpp2pRcS3YmNjSU9P97oMz2kql4iIzynoRUR8TkEvIuJznl2MNbPdwGZP3hwSOOKuXZ+Ltu8X9D1Hi2j8ngc55zodzxd4djHWOdfDq/c2s+zjvWodyaLt+wV9z9EiWr/n4/0aLd2IiPicgl5ExOeiNehnel1AC4u27xf0PUcLfc9N4NnFWBERaRnRekYvIhI1FPQiIj4XVUFvZpPMbIOZ5ZmZ7ydsmlmqmX1gZuvNbK2Z3eV1TS3FzGLMbKWZLfK6lpZgZl3N7FUz+7Lmv/eZXtfU3Mzs7pqf6zVmNtfMfLfhs5nNMrNCM1tT73Pdzew9M8ut+bNbY68TNUFfb5Pzi4AM4Dozy/C2qmZXCdzrnBsCnAHcEQXfc627qB6rHS1+B7ztnBsMDMPn37uZJQMzgEzn3FAghuq9MvzmBWDSEZ+7D3jfOTcAeJ8mjIWPmqCn3ibnzrlyoHaTc99yzu1wzq2o+Xgf1b/8yd5W1fzMLAW4GHjW61pagpl1Br5F9b4QOOfKnXN7vK2qRbQG2tXsateef975LuI55z7in3frmwzMrvl4NnB5Y68TTUHf0Cbnvg+9WmaWBowAlnpbSYv4LfBjIOh1IS2kL7AbeL5muepZM+vgdVHNyTm3DXgc2ALsAEqdc+96W1WLSXLO7YDqkzkgsbEviKagb9IG5n5kZh2B14B/dc7t9bqe5mRmlwCFzrnlXtfSgloDI4GnnXMjgAP4fJe3mnXpyUA60BvoYGbf9baq8BVNQd+UTc59x8xiqQ75Oc65v3pdTws4C7jMzDZRvTx3npm95G1Jza4AKHDO1f5t7VWqg9/Pzgc2Oud2O+cqgL8C4zyuqaXsMrNeADV/Fjb2BdEU9E3Z5NxXzMyoXrdd75z7jdf1tATn3H8651Kcc2lU/zf+X+ecr8/0nHM7ga1mNqjmUxOAdR6W1BK2AGeYWfuan/MJ+PwCdD0LgRtrPr4RWNDYF0TNVoLOuUozq93kPAaY5Zxb63FZze0sYCrwhZmtqvncT5xziz2sSZrHj4A5NScx+cDNHtfTrJxzS83sVWAF1d1lK/HhOAQzmwuMBxLMrAB4EHgUeMXMplH9P7wpjb6ORiCIiPhbNC3diIhEJQW9iIjPKehFRHxOQS8i4nMKehERn1PQi9Qws7T6UwJF/EJBLyLicwp6kQaYWd+aAWGjva5F5GQp6EWOUDNK4DXgZudcltf1iJysqBmBINJEPaieHXJlFIzIkCihM3qRbyqlet+Cs7wuRCRUdEYv8k3lVO/Y846Z7XfOvex1QSInS0EvcgTn3IGaDUzeM7MDzrlGx8CKhDNNrxQR8Tmt0YuI+JyCXkTE5xT0IiI+p6AXEfE5Bb2IiM8p6EVEfE5BLyLic/8feKtLKYGBYIsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "data.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We effectively have the following equations (I repeat the equation here for convenience):\n",
    "\n",
    "$$y(k) = a_1 y(k - 1) + a_2 y(k - 2) + b_1 u(k - 1) + b_2 u(k - 2)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.058 = a1×0.0 + a2×0.0 + b1×1 + b2×0\n",
      "0.22 = a1×0.058 + a2×0.0 + b1×1 + b2×1\n",
      "0.36 = a1×0.22 + a2×0.058 + b1×1 + b2×1\n",
      "0.49 = a1×0.36 + a2×0.22 + b1×1 + b2×1\n",
      "0.6 = a1×0.49 + a2×0.36 + b1×1 + b2×1\n",
      "0.69 = a1×0.6 + a2×0.49 + b1×1 + b2×1\n",
      "0.77 = a1×0.69 + a2×0.6 + b1×1 + b2×1\n",
      "0.83 = a1×0.77 + a2×0.69 + b1×1 + b2×1\n",
      "0.89 = a1×0.83 + a2×0.77 + b1×1 + b2×1\n",
      "0.93 = a1×0.89 + a2×0.83 + b1×1 + b2×1\n"
     ]
    }
   ],
   "source": [
    "for k in range(1, 11):\n",
    "    print(f'{y[k]:.2} = a1×{y[k - 1]:.2} + a2×{y[k - 2]:.2} + b1×{u[k - 1]} + b2×{u[k - 2]}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>k</th>\n",
       "      <th>y[k]</th>\n",
       "      <th>y[k-1]</th>\n",
       "      <th>y[k-2]</th>\n",
       "      <th>u[k-1]</th>\n",
       "      <th>u[k-2]</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>0.058</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>0.217</td>\n",
       "      <td>0.058</td>\n",
       "      <td>0.000</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>0.360</td>\n",
       "      <td>0.217</td>\n",
       "      <td>0.058</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4</td>\n",
       "      <td>0.488</td>\n",
       "      <td>0.360</td>\n",
       "      <td>0.217</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5</td>\n",
       "      <td>0.600</td>\n",
       "      <td>0.488</td>\n",
       "      <td>0.360</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>6</td>\n",
       "      <td>0.692</td>\n",
       "      <td>0.600</td>\n",
       "      <td>0.488</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>7</td>\n",
       "      <td>0.772</td>\n",
       "      <td>0.692</td>\n",
       "      <td>0.600</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>8</td>\n",
       "      <td>0.833</td>\n",
       "      <td>0.772</td>\n",
       "      <td>0.692</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>9</td>\n",
       "      <td>0.888</td>\n",
       "      <td>0.833</td>\n",
       "      <td>0.772</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>10</td>\n",
       "      <td>0.925</td>\n",
       "      <td>0.888</td>\n",
       "      <td>0.833</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    k   y[k]  y[k-1]  y[k-2]  u[k-1]  u[k-2]\n",
       "0   1  0.058   0.000   0.000       1       0\n",
       "1   2  0.217   0.058   0.000       1       1\n",
       "2   3  0.360   0.217   0.058       1       1\n",
       "3   4  0.488   0.360   0.217       1       1\n",
       "4   5  0.600   0.488   0.360       1       1\n",
       "5   6  0.692   0.600   0.488       1       1\n",
       "6   7  0.772   0.692   0.600       1       1\n",
       "7   8  0.833   0.772   0.692       1       1\n",
       "8   9  0.888   0.833   0.772       1       1\n",
       "9  10  0.925   0.888   0.833       1       1"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pandas.DataFrame([(k, \n",
    "                  y[k], \n",
    "                  y[k-1], \n",
    "                  y[k-2], \n",
    "                  u[k-1], \n",
    "                  u[k-2]) for k in range(1, 11)], columns=['k', 'y[k]', 'y[k-1]', 'y[k-2]', 'u[k-1]', 'u[k-2]'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We notice that these equations are linear in the coefficients. We can define $\\beta= [a_1, a_2, b_1, b_2]^T$. Now, to write the above equations in matrix form $Y = X \\beta $, we define"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.058],\n",
       "       [0.217],\n",
       "       [0.36 ],\n",
       "       [0.488],\n",
       "       [0.6  ],\n",
       "       [0.692],\n",
       "       [0.772],\n",
       "       [0.833],\n",
       "       [0.888],\n",
       "       [0.925]])"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Y = numpy.atleast_2d(y.loc[1:]).T\n",
    "Y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To build the coefficient matrix we observe that there are two blocks of constant diagonals (the part with the $y$s and the part with the $u$s). Matrices with constant diagonals are called Toeplitz matrices and can be constructed with the `scipy.linalg.toeplitz` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy\n",
    "import scipy.linalg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.   , 0.   ],\n",
       "       [0.058, 0.   ],\n",
       "       [0.217, 0.058],\n",
       "       [0.36 , 0.217],\n",
       "       [0.488, 0.36 ],\n",
       "       [0.6  , 0.488],\n",
       "       [0.692, 0.6  ],\n",
       "       [0.772, 0.692],\n",
       "       [0.833, 0.772],\n",
       "       [0.888, 0.833]])"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X1 = scipy.linalg.toeplitz(y.loc[0:9], [0, 0])\n",
    "X1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1, 0],\n",
       "       [1, 1],\n",
       "       [1, 1],\n",
       "       [1, 1],\n",
       "       [1, 1],\n",
       "       [1, 1],\n",
       "       [1, 1],\n",
       "       [1, 1],\n",
       "       [1, 1],\n",
       "       [1, 1]])"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X2 = scipy.linalg.toeplitz(u.loc[0:9], [0, 0])\n",
    "X2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = numpy.hstack([X1, X2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.   , 0.   , 1.   , 0.   ],\n",
       "       [0.058, 0.   , 1.   , 1.   ],\n",
       "       [0.217, 0.058, 1.   , 1.   ],\n",
       "       [0.36 , 0.217, 1.   , 1.   ],\n",
       "       [0.488, 0.36 , 1.   , 1.   ],\n",
       "       [0.6  , 0.488, 1.   , 1.   ],\n",
       "       [0.692, 0.6  , 1.   , 1.   ],\n",
       "       [0.772, 0.692, 1.   , 1.   ],\n",
       "       [0.833, 0.772, 1.   , 1.   ],\n",
       "       [0.888, 0.833, 1.   , 1.   ]])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another option is to use the loop from before to construct the matrices. This is a little more legible but slower:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y = []\n",
    "X = []\n",
    "for k in range(1, 11):\n",
    "    Y.append([y[k]])\n",
    "    X.append([y[k - 1], y[k - 2], u[k - 1], u[k - 2]])\n",
    "Y = numpy.array(Y)\n",
    "X = numpy.array(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We solve for $\\beta$ as we did for linear regression:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.98464753],\n",
       "       [-0.12211256],\n",
       "       [ 0.058     ],\n",
       "       [ 0.10124916]])"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "beta, _, _, _ = numpy.linalg.lstsq(X, Y, rcond=None)\n",
    "beta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}